# Load libraries
library(tidyverse)
library(tidylog)
library(sdmTMB)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(modelr)
library(ggstats)
library(ggsidekick); theme_set(theme_sleek())
home <- here::here()
source(paste0(home, "/R/functions/map-plot.R"))Fit example models
Read data and prediction grid, scales variables
# Read data
d <- readr::read_csv(paste0(home, "/data/clean/larval_size.csv")) %>%
drop_na(temp) %>%
mutate(temp_sc = (temp - mean(temp, na.rm = TRUE)) / sd(temp, na.rm = TRUE),
yday_ct = yday - mean(yday),
year_f = as.factor(year),
species_f = as.factor(species),
year_ct = year - median(year))
# Trim some outliers
d %>%
group_by(species) %>%
mutate(sd = sd(length_mm),
mean = mean(length_mm)) %>%
ungroup() %>%
mutate(outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
"No", "Yes")) %>%
ggplot(aes(year, length_mm, color = outlier)) +
facet_wrap(~species, scales = "free") +
geom_point()d %>%
group_by(species) %>%
mutate(sd = sd(length_mm),
mean = mean(length_mm)) %>%
ungroup() %>%
mutate(outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
"No", "Yes")) %>%
ggplot(aes(length_mm, fill = outlier)) +
facet_wrap(~species, scales = "free") +
geom_bar()sort(unique(d$species)) [1] "Agonus cataphractus" "Ammodytidae"
[3] "Anguilla anguilla" "Aphia minuta"
[5] "Argentina silus" "Argentina spp"
[7] "Chirolophis ascanii" "Clupea harengus"
[9] "Crystallogobius linearis" "Enchelyopus cimbrius"
[11] "Limanda limanda" "Maurolicus muelleri"
[13] "Microstomus kitt" "Myoxocephalus scorpius"
[15] "Pholis gunnellus" "Pomatoschistus sp"
[17] "Sardina pilchardus" "Sprattus sprattus"
[19] "Syngnathus rostellatus" "Taurulus bubalis"
d <- d %>%
group_by(species) %>%
mutate(sd = sd(length_mm),
mean = mean(length_mm)) %>%
ungroup() %>%
mutate(
outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
"No", "Yes"),
outlier = ifelse(species == "Limanda limanda" & length_mm > 50, "Yes", outlier)) %>%
filter(outlier == "No")
ggplot(d, aes(length_mm, fill = outlier)) +
facet_wrap(~species, scales = "free") +
geom_bar()d %>%
group_by(year, species) %>%
summarise(n = n()) %>%
ggplot(aes(year, n)) +
facet_wrap(~species, scales = "free_y", ncol = 3) +
geom_bar(stat = "identity")coul <- brewer.pal(11, "Spectral")
coul <- colorRampPalette(coul)(length(unique(d$species)))
p1 <- ggplot(d, aes(log(length_mm), fill = species)) +
geom_histogram() +
scale_fill_manual(values = coul, name = "Species") +
coord_cartesian(expand = 0) +
labs(y = "Count", x = "Length (mm)", tag = "b)") +
theme(legend.text = element_text(face = "italic", size = 7),
legend.key.size = unit(0.25, "cm"),
legend.position.inside = c(0.2, 0.69),
plot.tag = element_text()) +
guides(fill = guide_legend(position = "inside"))
# Load prediction grid
pred_grid <- readr::read_csv(paste0(home, "/data/clean/pred_grid.csv")) %>%
drop_na(temp) %>%
filter(year %in% unique(d$year)) %>%
mutate(temp_sc = (temp - mean(d$temp, na.rm = TRUE)) / sd(d$temp, na.rm = TRUE),
year_f = as.factor(year),
year_ct = 0,
yday_ct = 0) %>%
mutate(keep = ifelse(lon < 10 & lat < 57.15, "N", "Y")) %>%
filter(keep == "Y") %>% dplyr::select(-keep)
# Plot temperature in space
p2 <- plot_map +
geom_raster(data = pred_grid %>%
group_by(X, Y) %>%
summarise(mean_temp = mean(temp)),
aes(X*1000, Y*1000, fill = mean_temp)) +
geom_point(data = d %>% distinct(X, Y), aes(X*1000, Y*1000),
size = 0.5, alpha = 0.3, shape = 21, fill = NA, color = "white") +
scale_fill_viridis(option = "magma", name = "Temperature (°C)") +
theme(legend.position.inside = c(0.25, 0.1),
legend.direction = "horizontal",
legend.key.width = unit(0.5, "cm"),
legend.key.height = unit(0.2, "cm"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 8),
plot.tag = element_text()) +
labs(tag = "a)") +
guides(fill = guide_colorbar(position = "inside",
title.position = "top",
title.hjust = 0.5)) +
geom_sf() +
annotate("text", label = "Sweden", x = xmin2 + 0.95*xrange, y = ymin2 + 0.75*yrange,
color = "gray50", size = 3) +
annotate("text", label = "Norway", x = xmin2 + 0.08*xrange, y = ymin2 + 0.95*yrange,
color = "gray50", size = 3) +
annotate("text", label = "Denmark", x = xmin2 + 0.42*xrange, y = ymin2 + 0.45*yrange,
color = "gray50", size = 3)
p2 + p1ggsave(paste0(home, "/figures/data_map.pdf"), width = 22, height = 12, units = "cm")
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/00-fit-models_cache/html")Fit alternative models
mesh <- make_mesh(d,
xy_cols = c("X", "Y"),
cutoff = 5)
ggplot() +
inlabru::gg(mesh$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = d, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) +
labs(x = "Easting (km)", y = "Northing (km)")# Basic model
m0 <- sdmTMB(
length_mm ~ temp_sc*species + year_ct*species + yday_ct,
data = d,
mesh = mesh,
family = gengamma(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year")
sanity(m0)
m0Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc * species + year_ct * species + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: d
Family: gengamma(link = 'log')
coef.est coef.se
(Intercept) 2.18 0.02
temp_sc 0.03 0.01
speciesAmmodytidae 1.69 0.02
speciesAnguilla anguilla 2.09 0.02
speciesAphia minuta 1.34 0.01
speciesArgentina silus 1.58 0.16
speciesArgentina spp 1.28 0.03
speciesChirolophis ascanii 0.60 0.02
speciesClupea harengus 1.10 0.01
speciesCrystallogobius linearis 0.97 0.01
speciesEnchelyopus cimbrius 1.04 0.05
speciesLimanda limanda -0.39 0.03
speciesMaurolicus muelleri 1.46 0.04
speciesMicrostomus kitt 1.05 0.02
speciesMyoxocephalus scorpius -0.09 0.02
speciesPholis gunnellus 0.43 0.01
speciesPomatoschistus sp 1.49 0.01
speciesSardina pilchardus 1.30 0.08
speciesSprattus sprattus 1.72 0.02
speciesSyngnathus rostellatus 2.03 0.02
speciesTaurulus bubalis -0.22 0.05
year_ct 0.00 0.00
yday_ct 0.00 0.00
temp_sc:speciesAmmodytidae 0.02 0.01
temp_sc:speciesAnguilla anguilla -0.01 0.02
temp_sc:speciesAphia minuta -0.02 0.01
temp_sc:speciesArgentina silus -0.04 0.04
temp_sc:speciesArgentina spp -0.18 0.03
temp_sc:speciesChirolophis ascanii -0.02 0.01
temp_sc:speciesClupea harengus -0.05 0.01
temp_sc:speciesCrystallogobius linearis -0.03 0.01
temp_sc:speciesEnchelyopus cimbrius 0.19 0.05
temp_sc:speciesLimanda limanda 0.07 0.02
temp_sc:speciesMaurolicus muelleri -0.09 0.03
temp_sc:speciesMicrostomus kitt -0.11 0.02
temp_sc:speciesMyoxocephalus scorpius -0.04 0.02
temp_sc:speciesPholis gunnellus -0.06 0.01
temp_sc:speciesPomatoschistus sp -0.01 0.01
temp_sc:speciesSardina pilchardus -0.10 0.02
temp_sc:speciesSprattus sprattus -0.03 0.02
temp_sc:speciesSyngnathus rostellatus 0.00 0.01
temp_sc:speciesTaurulus bubalis -0.15 0.04
speciesAmmodytidae:year_ct -0.02 0.00
speciesAnguilla anguilla:year_ct 0.01 0.00
speciesAphia minuta:year_ct 0.00 0.00
speciesArgentina silus:year_ct -0.04 0.02
speciesArgentina spp:year_ct -0.01 0.00
speciesChirolophis ascanii:year_ct 0.00 0.00
speciesClupea harengus:year_ct 0.00 0.00
speciesCrystallogobius linearis:year_ct 0.00 0.00
speciesEnchelyopus cimbrius:year_ct -0.01 0.00
speciesLimanda limanda:year_ct 0.04 0.00
speciesMaurolicus muelleri:year_ct 0.02 0.00
speciesMicrostomus kitt:year_ct 0.00 0.00
speciesMyoxocephalus scorpius:year_ct 0.00 0.00
speciesPholis gunnellus:year_ct 0.00 0.00
speciesPomatoschistus sp:year_ct 0.01 0.00
speciesSardina pilchardus:year_ct -0.01 0.01
speciesSprattus sprattus:year_ct 0.00 0.00
speciesSyngnathus rostellatus:year_ct 0.00 0.00
speciesTaurulus bubalis:year_ct 0.01 0.01
Dispersion parameter: 0.22
Generalized gamma lambda: 0.51
Matérn range: 13.42
Spatial SD: 0.09
ML criterion at convergence: 149690.315
See ?tidy.sdmTMB to extract these values as a data frame.
# IID ST fields
m1 <- sdmTMB(
length_mm ~ temp_sc*species + year_ct*species + yday_ct,
data = d,
mesh = mesh,
family = gengamma(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year")
sanity(m1)
m1Spatiotemporal model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc * species + year_ct * species + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: d
Family: gengamma(link = 'log')
coef.est coef.se
(Intercept) 2.20 0.02
temp_sc 0.02 0.01
speciesAmmodytidae 1.65 0.02
speciesAnguilla anguilla 2.08 0.02
speciesAphia minuta 1.31 0.01
speciesArgentina silus 1.64 0.17
speciesArgentina spp 1.27 0.03
speciesChirolophis ascanii 0.58 0.02
speciesClupea harengus 1.09 0.01
speciesCrystallogobius linearis 0.96 0.01
speciesEnchelyopus cimbrius 0.99 0.05
speciesLimanda limanda -0.34 0.03
speciesMaurolicus muelleri 1.43 0.04
speciesMicrostomus kitt 1.05 0.02
speciesMyoxocephalus scorpius -0.08 0.02
speciesPholis gunnellus 0.40 0.01
speciesPomatoschistus sp 1.47 0.01
speciesSardina pilchardus 1.26 0.09
speciesSprattus sprattus 1.70 0.02
speciesSyngnathus rostellatus 2.01 0.02
speciesTaurulus bubalis -0.21 0.05
year_ct -0.01 0.00
yday_ct 0.00 0.00
temp_sc:speciesAmmodytidae 0.03 0.01
temp_sc:speciesAnguilla anguilla -0.01 0.02
temp_sc:speciesAphia minuta -0.02 0.01
temp_sc:speciesArgentina silus -0.03 0.04
temp_sc:speciesArgentina spp -0.17 0.03
temp_sc:speciesChirolophis ascanii -0.02 0.01
temp_sc:speciesClupea harengus -0.05 0.01
temp_sc:speciesCrystallogobius linearis -0.03 0.01
temp_sc:speciesEnchelyopus cimbrius 0.17 0.05
temp_sc:speciesLimanda limanda 0.08 0.02
temp_sc:speciesMaurolicus muelleri -0.10 0.03
temp_sc:speciesMicrostomus kitt -0.12 0.02
temp_sc:speciesMyoxocephalus scorpius -0.03 0.02
temp_sc:speciesPholis gunnellus -0.05 0.01
temp_sc:speciesPomatoschistus sp -0.02 0.01
temp_sc:speciesSardina pilchardus -0.10 0.02
temp_sc:speciesSprattus sprattus -0.03 0.02
temp_sc:speciesSyngnathus rostellatus -0.01 0.01
temp_sc:speciesTaurulus bubalis -0.12 0.04
speciesAmmodytidae:year_ct -0.02 0.00
speciesAnguilla anguilla:year_ct 0.01 0.00
speciesAphia minuta:year_ct 0.00 0.00
speciesArgentina silus:year_ct -0.05 0.02
speciesArgentina spp:year_ct -0.01 0.00
speciesChirolophis ascanii:year_ct 0.00 0.00
speciesClupea harengus:year_ct 0.00 0.00
speciesCrystallogobius linearis:year_ct 0.00 0.00
speciesEnchelyopus cimbrius:year_ct 0.00 0.00
speciesLimanda limanda:year_ct 0.03 0.00
speciesMaurolicus muelleri:year_ct 0.02 0.00
speciesMicrostomus kitt:year_ct 0.00 0.00
speciesMyoxocephalus scorpius:year_ct 0.00 0.00
speciesPholis gunnellus:year_ct 0.01 0.00
speciesPomatoschistus sp:year_ct 0.01 0.00
speciesSardina pilchardus:year_ct -0.01 0.01
speciesSprattus sprattus:year_ct 0.00 0.00
speciesSyngnathus rostellatus:year_ct 0.00 0.00
speciesTaurulus bubalis:year_ct 0.00 0.01
Dispersion parameter: 0.21
Generalized gamma lambda: 0.51
Matérn range: 51.48
Spatial SD: 0.04
Spatiotemporal IID SD: 0.09
ML criterion at convergence: 148310.742
See ?tidy.sdmTMB to extract these values as a data frame.
Cross validation
set.seed(99)
clust <- kmeans(d[, c("X", "Y")], 5)$cluster
d$clust <- clust
pal <- coul[c(1, 6, 8, 16, 20)]
plot_map +
geom_point(data = d, aes(X*1000, Y*1000, color = as.factor(clust)),
alpha = 0.3, size = 0.6, shape = 21, fill = NA) +
scale_color_manual(values = pal)Warning: `guide_colourbar()` needs continuous scales.
ggsave(paste0(home, "/figures/cv_cluster.pdf"), width = 11, height = 11, units = "cm")Warning: `guide_colourbar()` needs continuous scales.
m0_cv <- sdmTMB_cv(
length_mm ~ temp_sc*species + year_ct*species + yday_ct,
data = d,
mesh = mesh,
family = gengamma(link = "log"),
spatiotemporal = "off",
spatial = "on",
fold_ids = clust,
k_folds = length(unique(clust))
)
m1_cv <- sdmTMB_cv(
length_mm ~ temp_sc*species + year_ct*species + yday_ct,
data = d,
mesh = mesh,
family = gengamma(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
fold_ids = clust,
k_folds = length(unique(clust))
)
# Compare log-likelihoods -- higher is better!
m0_cv$sum_loglik[1] -154327.9
m1_cv$sum_loglik[1] -153179.6
Go with model 1, spatiotemporal IID
Plot residuals
# samps <- sdmTMBextra::predict_mle_mcmc(m1, mcmc_iter = 201, mcmc_warmup = 200)
# d$mcmc_res1 <- as.vector(residuals(m1, type = "mle-mcmc", mcmc_samples = samps))
# #d$mcmc_res1 <- residuals(m1)
#
# # Looks reasonably, a weird at the tails though
# d %>%
# ggplot(aes(sample = mcmc_res1)) +
# stat_qq(size = 0.75, shape = 21, fill = NA) +
# stat_qq_line() +
# labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
# theme(aspect.ratio = 1)
#
# ggsave(paste0(home, "/figures/qq.pdf"), width = 11, height = 11, units = "cm")
#
# # Plot in space
# plot_map +
# geom_point(data = d, aes(X*1000, Y*1000, color = mcmc_res1)) +
# scale_color_gradient2()
# Residuals -- above doesn't work yet with gengamma!
png(file = paste0(home, "/figures/qq.png"), units = "cm", width = 15, height = 15, res = 300)
simulate(m1, nsim = 300, type = 'mle-mvn') %>%
dharma_residuals(m1)
dev.off()quartz_off_screen
2
Predict on grid and plot spatial
pred <- predict(m1, newdata = pred_grid %>%
mutate(species = "Clupea harengus")) # The ST and S random fields do not change with species
# Spatial random
plot_map +
geom_raster(data = pred %>% filter(year == max(year)),
aes(X*1000, Y*1000, fill = omega_s)) +
scale_fill_gradient2(name = "Spatial random effect") +
theme_sleek(base_size = 9) +
theme(legend.position.inside = c(0.25, 0.1),
legend.direction = "horizontal",
legend.key.width = unit(0.7, "cm"),
legend.key.height = unit(0.3, "cm")) +
guides(fill = guide_colorbar(position = "inside",
title.position = "top",
title.hjust = 0.5)) +
geom_sf()ggsave(paste0(home, "/figures/omega.pdf"), width = 11, height = 11, units = "cm")
# Spatiotemporal random
plot_map_fc +
geom_raster(data = pred,
aes(X*1000, Y*1000, fill = epsilon_st)) +
scale_fill_gradient2(name = "Spatiotemporal random effects") +
facet_wrap(~year, ncol = 7) +
theme_facet_map(base_size = 10) +
geom_sf()ggsave(paste0(home, "/figures/epsilon.pdf"), width = 17, height = 15, units = "cm")
# Prediction (herring), selected years
plot_map_fc +
geom_raster(data = pred %>% filter(year %in% c(1992, 2002, 2012, 2022)),
aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(name = "Length (mm)", trans = "sqrt") +
facet_wrap(~year, ncol = 4) +
theme_sleek(base_size = 9) +
theme(legend.position.inside = c(0.195, 0.75),
legend.direction = "vertical",
legend.key.width = unit(0.2, "cm"),
legend.key.height = unit(0.3, "cm"),
legend.title = element_text(size = 6)) +
guides(fill = guide_colorbar(position = "inside",
title.position = "top",
title.hjust = 1)) +
geom_sf()ggsave(paste0(home, "/figures/st_pred.pdf"), width = 17, height = 7, units = "cm")
# Plot uncertainty
sim <- predict(m1, nsim = 100,
newdata = pred_grid %>%
mutate(species = "Clupea harengus"))
sim_last <- sim[pred_grid$year == max(pred_grid$year), ] # just plot last year
pred_last <- pred[pred$year == max(pred_grid$year), ]
pred_last$lwr <- apply(exp(sim_last), 1, quantile, probs = 0.025)
pred_last$upr <- apply(exp(sim_last), 1, quantile, probs = 0.975)
pred_last$sd <- round(apply(exp(sim_last), 1, function(x) sd(x)), 2)
pred_last$cv <- round(apply(exp(sim_last), 1, function(x) sd(x) / mean(x)), 2)
plot_map +
geom_raster(data = pred_last,
aes(X*1000, Y*1000, fill = cv)) +
scale_fill_viridis(option = "mako") +
theme_sleek(base_size = 9) +
theme(legend.position.inside = c(0.25, 0.1),
legend.direction = "horizontal",
legend.key.width = unit(0.7, "cm"),
legend.key.height = unit(0.3, "cm")) +
guides(fill = guide_colorbar(position = "inside",
title.position = "top",
title.hjust = 0.5)) +
geom_sf()ggsave(paste0(home, "/figures/spatial_cv.pdf"), width = 11, height = 11, units = "cm")Coefficients
tidy(m1) %>% as.data.frame() term estimate std.error
1 (Intercept) 2.2017933324 0.0170624303
2 temp_sc 0.0246658086 0.0128297468
3 speciesAmmodytidae 1.6468346729 0.0160873369
4 speciesAnguilla anguilla 2.0773744737 0.0194810177
5 speciesAphia minuta 1.3087453316 0.0136375644
6 speciesArgentina silus 1.6408805276 0.1678850380
7 speciesArgentina spp 1.2661382522 0.0303506268
8 speciesChirolophis ascanii 0.5844242365 0.0157352162
9 speciesClupea harengus 1.0942305689 0.0135770762
10 speciesCrystallogobius linearis 0.9636805608 0.0136184992
11 speciesEnchelyopus cimbrius 0.9921311106 0.0469485503
12 speciesLimanda limanda -0.3413526086 0.0307652385
13 speciesMaurolicus muelleri 1.4316224155 0.0400275537
14 speciesMicrostomus kitt 1.0468385677 0.0165714706
15 speciesMyoxocephalus scorpius -0.0795154997 0.0191302853
16 speciesPholis gunnellus 0.4031977450 0.0136823916
17 speciesPomatoschistus sp 1.4651207322 0.0142411280
18 speciesSardina pilchardus 1.2560352585 0.0909112746
19 speciesSprattus sprattus 1.6979418847 0.0176364136
20 speciesSyngnathus rostellatus 2.0109400867 0.0157665100
21 speciesTaurulus bubalis -0.2103747533 0.0521629182
22 year_ct -0.0051828622 0.0013064920
23 yday_ct 0.0018833563 0.0005274263
24 temp_sc:speciesAmmodytidae 0.0312250113 0.0138272361
25 temp_sc:speciesAnguilla anguilla -0.0072353640 0.0202931399
26 temp_sc:speciesAphia minuta -0.0204725003 0.0121542931
27 temp_sc:speciesArgentina silus -0.0258507957 0.0374055255
28 temp_sc:speciesArgentina spp -0.1744859525 0.0335999122
29 temp_sc:speciesChirolophis ascanii -0.0226885606 0.0136453105
30 temp_sc:speciesClupea harengus -0.0493219672 0.0123029235
31 temp_sc:speciesCrystallogobius linearis -0.0334516293 0.0121512676
32 temp_sc:speciesEnchelyopus cimbrius 0.1705620970 0.0459767541
33 temp_sc:speciesLimanda limanda 0.0800843872 0.0184944242
34 temp_sc:speciesMaurolicus muelleri -0.0958958240 0.0316982103
35 temp_sc:speciesMicrostomus kitt -0.1187109721 0.0168959849
36 temp_sc:speciesMyoxocephalus scorpius -0.0270361945 0.0152627963
37 temp_sc:speciesPholis gunnellus -0.0534576875 0.0120790401
38 temp_sc:speciesPomatoschistus sp -0.0167895195 0.0125776493
39 temp_sc:speciesSardina pilchardus -0.1020761875 0.0190106294
40 temp_sc:speciesSprattus sprattus -0.0323183759 0.0155200946
41 temp_sc:speciesSyngnathus rostellatus -0.0137953804 0.0137118139
42 temp_sc:speciesTaurulus bubalis -0.1201217321 0.0429687448
43 speciesAmmodytidae:year_ct -0.0169662389 0.0014053404
44 speciesAnguilla anguilla:year_ct 0.0053383930 0.0019758344
45 speciesAphia minuta:year_ct 0.0015989702 0.0012125807
46 speciesArgentina silus:year_ct -0.0497146976 0.0176120858
47 speciesArgentina spp:year_ct -0.0079958455 0.0022274962
48 speciesChirolophis ascanii:year_ct -0.0012728807 0.0013582374
49 speciesClupea harengus:year_ct 0.0015120757 0.0011879410
50 speciesCrystallogobius linearis:year_ct -0.0002908017 0.0011894349
51 speciesEnchelyopus cimbrius:year_ct -0.0013727965 0.0034116001
52 speciesLimanda limanda:year_ct 0.0320642189 0.0027933830
53 speciesMaurolicus muelleri:year_ct 0.0246322658 0.0035843246
54 speciesMicrostomus kitt:year_ct -0.0009833155 0.0015857900
55 speciesMyoxocephalus scorpius:year_ct -0.0002173287 0.0018444513
56 speciesPholis gunnellus:year_ct 0.0054458029 0.0012005775
57 speciesPomatoschistus sp:year_ct 0.0077844893 0.0012831516
58 speciesSardina pilchardus:year_ct -0.0090907467 0.0059110581
59 speciesSprattus sprattus:year_ct -0.0024785831 0.0017975039
60 speciesSyngnathus rostellatus:year_ct 0.0006984562 0.0015935258
61 speciesTaurulus bubalis:year_ct 0.0042812992 0.0087080424
sort(unique(d$species)) [1] "Agonus cataphractus" "Ammodytidae"
[3] "Anguilla anguilla" "Aphia minuta"
[5] "Argentina silus" "Argentina spp"
[7] "Chirolophis ascanii" "Clupea harengus"
[9] "Crystallogobius linearis" "Enchelyopus cimbrius"
[11] "Limanda limanda" "Maurolicus muelleri"
[13] "Microstomus kitt" "Myoxocephalus scorpius"
[15] "Pholis gunnellus" "Pomatoschistus sp"
[17] "Sardina pilchardus" "Sprattus sprattus"
[19] "Syngnathus rostellatus" "Taurulus bubalis"
mean_coefs <- tidy(m1, conf.int = TRUE) %>%
filter(term %in% c("temp_sc", "yday_ct")) %>%
mutate(term = ifelse(term == "temp_sc",
"temp_sc:speciesAgonus cataphractus",
"speciesAgonus cataphractus:year_ct"))
# Wrangle data
coefs_temp <- tidy(m1, conf.int = TRUE) %>%
filter(grepl("temp_sc:", term)) %>%
bind_rows(mean_coefs %>% filter(term == "temp_sc:speciesAgonus cataphractus")) %>%
mutate(estimate2 = mean_coefs %>%
filter(term == "temp_sc:speciesAgonus cataphractus") %>%
pull(estimate) + estimate,
conf.low2 = mean_coefs %>%
filter(term == "temp_sc:speciesAgonus cataphractus") %>%
pull(conf.low) + conf.low,
conf.high2 = mean_coefs %>%
filter(term == "temp_sc:speciesAgonus cataphractus") %>%
pull(conf.high) + conf.high) %>%
mutate(species = str_remove(term, "temp_sc:"),
species = str_remove(species, "species"))
coefs_year <- tidy(m1, conf.int = TRUE) %>%
filter(grepl(":year_ct", term)) %>%
bind_rows(mean_coefs %>% filter(term == "speciesAgonus cataphractus:year_ct")) %>%
mutate(estimate2 = mean_coefs %>%
filter(term == "speciesAgonus cataphractus:year_ct") %>%
pull(estimate) + estimate,
conf.low2 = mean_coefs %>%
filter(term == "speciesAgonus cataphractus:year_ct") %>%
pull(conf.low) + conf.low,
conf.high2 = mean_coefs %>%
filter(term == "speciesAgonus cataphractus:year_ct") %>%
pull(conf.high) + conf.high) %>%
mutate(species = str_remove(term, ":year_ct"),
species = str_remove(species, "species"))
p1 <- coefs_temp %>%
mutate(sign = ifelse(estimate2 > 0, "pos", "neg"),
sig = ifelse(estimate2 > 0 & conf.low2 > 0, "sig", "not sig"),
sig = ifelse(estimate2 < 0 & conf.high2 < 0, "sig", sig)) %>%
ggplot(aes(estimate2, reorder(species, coefs_year$estimate2), fill = sig, color = sign, shape = sig)) +
geom_point(fill = NA) +
scale_shape_manual(values = c(21, 19)) +
geom_errorbar(aes(xmin = conf.low2, xmax = conf.high2),
width = 0, alpha = 0.3) +
geom_vline(xintercept = 0, alpha = 0.3, linetype = 2, linewidth = 0.25) +
theme(axis.text.y = element_text(face = "italic")) +
labs(y = "Species", x = "Temperature slope") +
scale_color_brewer(palette = "Dark2", direction = -1) +
scale_fill_brewer(palette = "Dark2", direction = -1) +
guides(color = "none", shape = "none") +
geom_stripped_rows(aes(y = species), inherit.aes = FALSE) +
theme(legend.position = "bottom")
p2 <- coefs_year %>%
mutate(sign = ifelse(estimate2 > 0, "pos", "neg"),
sig = ifelse(estimate2 > 0 & conf.low2 > 0, "sig", "not sig"),
sig = ifelse(estimate2 < 0 & conf.high2 < 0, "sig", sig)) %>%
ggplot(aes(estimate2, reorder(species, coefs_year$estimate2), fill = sig, color = sign, shape = sig)) +
geom_point(fill = NA) +
scale_shape_manual(values = c(21, 19)) +
geom_errorbar(aes(xmin = conf.low2, xmax = conf.high2),
width = 0, alpha = 0.3) +
geom_vline(xintercept = 0, alpha = 0.3, linetype = 2, linewidth = 0.25) +
theme(axis.text.y = element_text(face = "italic")) +
labs(y = "Species", x = "Year slope") +
scale_color_brewer(palette = "Dark2", direction = -1) +
scale_fill_brewer(palette = "Dark2", direction = -1) +
guides(color = "none", shape = "none") +
geom_stripped_rows(aes(y = species), inherit.aes = FALSE) +
theme(legend.position = "bottom")
p2 + p1 + plot_layout(axes = "collect")ggsave(paste0(home, "/figures/coefs.pdf"), width = 17, height = 10, units = "cm")Conditional predictions of temperature and year, with real values on the x and data, on 3 select species with contrasting trends
# Year effect
nd_y <- d %>%
filter(species %in% c("Limanda limanda", "Maurolicus muelleri", "Sardina pilchardus")) %>%
group_by(species) %>%
data_grid(year_ct = unique(year_ct)) %>%
ungroup() %>%
mutate(yday_ct = 0,
temp_sc = 0,
year = as.numeric(year_ct + median(d$year))) %>%
filter(year %in% unique(d$year))
pred_y <- predict(m1, newdata = nd_y, re_form = NA, re_form_iid = NA, se_fit = TRUE)
# Temperature effect
nd_t <- d %>%
filter(species %in% c("Limanda limanda", "Maurolicus muelleri", "Sardina pilchardus")) %>%
group_by(species) %>%
data_grid(temp_sc = unique(temp_sc)) %>%
ungroup() %>%
mutate(yday_ct = 0,
year_ct = 10, # 2015
year = as.numeric(year_ct + median(d$year)),
temp = (temp_sc * sd(d$temp)) + mean(d$temp)) %>%
filter(year %in% unique(d$year))
pred_t <- predict(m1, newdata = nd_t, re_form = NA, re_form_iid = NA, se_fit = TRUE)
# Plot!
lab_y <- tibble(species = c("Limanda limanda", "Maurolicus muelleri", "Sardina pilchardus"),
lab = c("a)", "b)", "c)"))
p1 <- ggplot(pred_y, aes(year, exp(est), color = species, fill = species)) +
geom_jitter(data = d %>% filter(species %in% c(nd_y$species)),
aes(year, length_mm), alpha = 0.5, height = 0, width = 0.3,
size = 0.95, fill = NA, shape = 21) +
geom_ribbon(aes(ymin = exp(est - 1.96*est_se),
ymax = exp(est + 1.96*est_se)), alpha = 0.4,
fill = "grey", color = NA) +
geom_line(linewidth = 1, color = "grey50") +
theme(strip.text = element_text(face = "italic")) +
geom_text(data = lab_y, aes(label = lab), x = 1993, y = Inf, vjust = 2,
color = "grey20") +
labs(x = "Year", y = "Length (mm)") +
coord_cartesian(xlim = c(min(pred_y$year), max(pred_y$year))) +
guides(color = "none", fill = "none") +
facet_wrap(~species, scales = "free") +
scale_color_brewer(palette = "Accent") +
scale_fill_brewer(palette = "Accent")
p1lab_t <- tibble(species = c("Limanda limanda", "Maurolicus muelleri", "Sardina pilchardus"),
lab = c("d)", "e)", "f)"))
p2 <- ggplot(pred_t, aes(temp, exp(est), color = species, fill = species)) +
geom_jitter(data = d %>% filter(species %in% c(nd_y$species)),
aes(temp, length_mm), alpha = 0.5, height = 0, width = 0.3,
size = 0.95, fill = NA, shape = 21) +
geom_ribbon(aes(ymin = exp(est - 1.96*est_se),
ymax = exp(est + 1.96*est_se)), alpha = 0.4,
fill = "grey", color = NA) +
geom_line(linewidth = 1, color = "grey50") +
theme(strip.text = element_text(face = "italic")) +
geom_text(data = lab_t, aes(label = lab), x = 0.5, y = Inf, vjust = 2,
color = "grey20") +
labs(x = "Temperature (°C)", y = "Length (mm)") +
coord_cartesian(xlim = c(min(pred_t$temp), max(pred_t$temp))) +
guides(color = "none", fill = "none") +
facet_wrap(~species, scales = "free") +
scale_color_brewer(palette = "Accent") +
scale_fill_brewer(palette = "Accent")
p2p1 / p2 +
plot_layout(axes = "collect")ggsave(paste0(home, "/figures/conditional.pdf"), width = 17, height = 12, units = "cm")